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Abstract. In this work we calculate the thermodynamic properties of hydrogen- 
helium plasmas with different mass fractions of helium by the direct path integral 
Monte Carlo method. To avoid unphysical approximations we use the path integral 
representation of the density matrix. We pay special attention to the region of weak 
coupling and degeneracy and compare the results of simulation with a model based 
on the chemical picture. Further with the help of calculated deuterium isochors we 
compute the shock Hugoniot of deuterium. We analyze our results in comparison with 
recent experimental and calculated data on the deuterium Hugoniot. 
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1. Introduction 

Hydrogen and helium are the most abundant elements in the Universe, therefore 
thermodynamic properties of hydrogen and helium plasmas are widely required for 
many astrophysical problems [1-4]. In particular, the investigation of the giant planets 
Jupiter and Saturn, and to a lesser extent brown dwarfs demands the thermodynamic 
information for hydrogen and helium in the approximate range of temperatures 
10^ < T < 10^ K and mass densities 0.01 < p < 100 g/cm^. This region is 
characterized by coupling effects and chemical reactions caused by partial pressure 
dissociation and ionization [5, 6]; these effects considerably complicate an equation 
of state (EOS) calculation. Moreover, in the same range of parameters the so-called 
plasma phase transition (PPT) has been predicted by many authors [1,5-8]. However 
the apphcation of the chemical picture [5, 6] at densities corresponding to pressure 
ionization is questionable. Therefore there is a great interest in direct first-principle 
numerical simulations of strongly coupled degenerate systems which avoid difficulties of 
conventional theories. 
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In this work we use the direct path integral Monte Carlo method (DPIMC) to 
calculate thermodynamic properties of hydrogen-helium plasma with different mass 
fractions of helium. This method is well established theoretically and allows the 
treatment of quantum and exchange effects without any approximations using only 
fundamental physical constants. We compare the results of our simulation with the 
EOS model based on the chemical picture [1,2]. We also use the DPIMC method 
to compute the deuterium Hugoniot. We compare our simulation results with recent 
experimental and theoretical works and analyze the modern state of the problem. 

2. Simulation method and results for hydrogen-helium plasma 

The details of our computational scheme can be found elsewhere [9-12]. Modern 
supercomputers allow us to simulate about 100 quantum particles in a Monte Carlo 
cell at a given temperature and volume. The DPIMC has no limitations on coupling 
parameter and can be applied at significant degeneracy of the system (with degeneracy 
parameter values as high as 300) [10]. Earlier the method was thoroughly tested by 
simulating different properties of ideal and interacting degenerate plasmas [13,14]. In 
particular, we investigated temperature and pressure dissociation and ionization ab 
initio; we also observed the effect of proton ordering at very high densities and the 
formation of a Coulomb crystal of protons [13]. 

In this section we calculate thermodynamic properties of hydrogen-helium mixtures 
at relatively low coupling and degeneracy parameters and compare our results with 
a well-known chemical picture model used mostly in astrophysics [1,2]. This model 
includes classical statistics for molecules and ions and Fermi-Dirac statistics for the 
electrons. It takes into account many physical effects including a number of subtle 
"second-order" phenomena. We calculated thermodynamic properties of hydrogen- 
helium mixtures with a composition corresponding to that of the outer layers of the 
Jovian atmosphere. During the mission of the Galileo spacecraft the helium abundance 
in the atmosphere of Jupiter was determined as F = mne/ (^He + '^h) = 0.234 and was 
close to the present-day protosolar value Y = 0.275 [3]. As the model of the Jupiter is 
significantly determined by its composition and EOS, it was interesting to simulate the 
thermodynamic properties of the mixture with different compositions in the region of 
pressure dissociation and ionization. 

We considered two mixtures with low and high abundance of helium. The results of 
calculations for the mixture corresponding to the outer layers of the Jovian atmosphere 
{Y = 0.234) in the region of temperatures from T = 10^ to 2 • 10^ K and electron number 
densities from Ue = 10^^ to 3-10^^ cm~^ are presented in figure^ The agreement between 
our calculations and the model [2] along the isotherms T = 4-10^, 5-10^, 10^, and 2-10^ K 
is quite good and becomes better with the increase of temperature. The formation of 
atoms and molecules is the reason of the pressure and energy reduction along the 10^ K 
isotherm with respect to the isotherm of a non-interacting hydrogen-helium mixture 
(see figure Hj). 
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Figure 1. Pressure (a) and energy per particle (b) in a liydrogen-licliuni mixture 
with the mass concentration of helium Y = 0.234 (Ry « 13.6 eV). Shown arc DPIMC 
isotherms and related EOS isotherms [2]. EOS [2] (DPIMC) calculations: 1(5) - 40 kK, 
2(6) — 50 kK, 3(7) — 100 kK, 4(8) — 200 kK. 9 - 100 kK isotherm for ideal plasma. 




Figure 2. Pressure (a) and energy per particle (b) in a hydrogen-helium mixture 
with the mass concentration of helium Y = 0.988. Shown are DPIMC isotherms and 
related EOS isotherms [2]. EOS [2] (DPIMC) calculations: 1(5) — 100 kK, 2(6) — 
156 kK, 3(7) — 200 kK, 4(8) — 312 kK. 9 — 100 kK isotherm for ideal plasma. C2 — 
critical point of the PPT [8] {Tcr « 120 kK). 



The results for Y = 0.988 (almost pure helium) at relatively high temperatures 
T = 10^ -V- 3 • 10^ K in a wide range of densities are presented in figure |21 The agreement 
between our calculations and the model [2] along the isotherms T = 10^, 1.56-10^, 2-10^, 
and 3.12 • 10^ K is satisfactory for pressure and internal energy per particle. The smaller 
values of pressure on the DPIMC isotherms 10^ and 1.56 ■ 10^ K near the particle density 
10^^ cm~^ can be explained by a strong influence of interaction and bound states in this 
region (see below). Ionization effects also reduce the internal energy of the system in 
comparison with non-interacting (ideal) plasma as it can be clearly seen in Fig. |2b- The 
positions of ionization minima are well reproduced by the DPIMC method in a good 
agreement with the chemical picture calculations. At higher densities Fermi-repulsion 
gives the main contribution to pressure and energy and this effect is also observed in 
our simulations. 

At low temperatures T < 3 • 10^ K and Y = 0.234 the agreement between 
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DPIMC and chemical picture calculations becomes worse, moreover, the region of 
thermodynamic instability has been discovered. In particular, along the isotherm 
T = 2 ■ 10^ K we have found such a region in the range of densities between 0.5 and 
5 g/ cm^. Along the isotherms T = 1.5-10^ K and T = 10^ K this region is even wider and 
begins from 0.38 g/cm^ [12]. Surprisingly, the region of DPIMC instability correlates 
with the range of temperatures (T < 2 ■ 10^ K) and densities (0.3-1 g/cm'^) in which 
the PPT in hydrogen or hydrogen-helium mixture with low mass concentration has 
been predicted [1,8,15]. Moreover, the sharp rise of electrical conductivity of hydrogen- 
helium mixture along the quasi-isentrope is also revealed experimentally in the range 
of densities 0.5-0.83 g/cm^ [16]. However, we cannot claim that these facts confirm the 
existence of PPT in our DPIMC simulation; in the nearest future we plan to investigate 
the PPT problem in detail using more sophisticated numerical methods. 

Because of the high binding energy of electrons in He we currently can obtain 
reliable results for Y = 0.988 only at temperatures higher than 10^ K. Under these 
conditions the influence of helium double ionization can lead to the formation of bound 
states in the Monte Carlo cell as well as pressure and internal energy decrease. Probably 
this effect takes place in Fig. linear electron number density Ue = 10^^ cm~^ at T = 10^ 
and 1.56 • 10^ K; the critical point of the possible PPT in this region with critical 
temperature ~ 120 kK [8] is also shown in figure Eti- 

3. Deuterium shock Hugoniot 

Using our previous simulation results for deuterium we calculated the shock 
Hugoniot of liquid deuterium [17]. Figure El summarizes the data from different 
experimental, theoretical, and numerical studies on the shock compression of deuterium. 
Measurements performed in the NOVA facility, where a shock wave in liquid deuterium 
with initial density 0.171 g/cm^ was generated by a laser pulse [18, 19] show that the 
deuterium density behind the shock front can increase by a factor of more than 6. 
Experiments with the acceleration of an aluminum foil by a magnetic field to velocities 
higher than 20 km/s [20] show a considerably lower compression ratio in comparison 
to [18,19]. The results obtained in [18,19] and [20] disagree within experimental errors. 

In contrast to [18,19] and [20], where targets several hundred microns thick were 
used, in [21-24], the shock compressibility of solid (initial density 0.199 g/cm^) [21-23] 
and liquid [23] deuterium was measured in a 4-mm-thick layer using a hemispherical 
explosive device. It is interesting to note that first such measurements for solid 
deuterium [21,23] (points 4) showed greater compressibility of deuterium than it was 
reported later [24] (points 5 in figure El). The same situation is observed for the 
experimental points on liquid deuterium (points 6 and 7, correspondingly, see [25] where 
preliminary experimental data for liquid deuterium from [24] are shown). Experimental 
points for liquid deuterium [24] are in a good correspondence with the data [20]. Another 
hemispherical device was applied for shock loading of dense gaseous deuterium with 
initial density close to that of liquid deuterium [25] . In these experiments [25] apart from 
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Figures. Shock Hugoniot of deuterium. Experimental data for liquid deuterium: 1 — 
[18,19], 2 — [20], 3 — [26], 6 — [25], and 7 — [24]; for solid deuterium: 4 — [21,23], 
and 5 — [24]; for gaseous deuterium: 8 — [25]. Calculations: 9 — [27], 10 — [28], 
11 — [29], 12 — [30], 13 — [31], 14 — [32], 15 — [25], 16 — [33], and 17 — this study. 

kinematic shock wave parameters temperature and light absorption of shock-compressed 
gas were registered. Two experimental points 8 corresponding to the initial gas densities 
0.1335 g/cm^ and 0.153 g/cm^ are also shown in figure El Curve 15 demonstrates 
the SAHA-IV liquid deuterium Hugoniot with the initial density 0.171 g/cm^ [25]. 
The SAHA-IV chemical plasma model was calibrated so as to be in agreement with 
points 8. In this case curve 15 passes through the old position of the liquid Hugoniot 
point at 1.09 Mbar [25]. The new position of the point at 1.09 Mbar [24], however, is 
shifted towards lower densities. Therefore points 7 and 8 in figure El probably cannot be 
described by one and the same theoretical model. 

In figure El a number of calculated shock Hugoniots is also shown; the detailed 
analysis of these results can be found in our recent works [17,34]. Here we can only 
indicate that the DPIMC Hugoniot is shifted towards higher densities in comparison 
to the experimental data published in [20,24]. At pressures below 1-2 Mbar, the 
thermodynamic instability revealed in [35] comes into play; therefore, a segment of 
the shock Hugoniot that lies below 1 Mbar is not quite reliable. At higher pressures the 
closest to the DPIMC Hugoniot is curve 16 calculated in [33] by the classical reactive 
ensemble Monte Carlo method. In this method, the effects of dissociation of deuterium 
molecules are taken into account most correctly; this allows one to achieve good 
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agreement with the experimental data obtained at low temperatures and pressures [26], 
even if ionization is not taken into account. Therefore we combined the low-pressure 
part of the Hugoniot from [33] and high-pressure one from [17] at 15000 K and obtained 
the united Hugoniot [34] (curve 17 in figure EI). We want to stress here that these two 
methods are completely independent and no interpolation procedure is used. 

Thus we confirm that the experimental points [18, 19] are questionable and the 
true position of the liquid deuterium Hugoniot remains unclear. We believe that 
future experiments at the hemispherical device [25] for densities of gaseous deuterium 
corresponding to the liquid and solid states will give important additional information 
about shock compression of liquid and solid deuterium. In the nearest future we plan to 
calculate two DPIMC Hugoniots corresponding to the initial gaseous deuterium densities 
0.1335 and 0.153 g/cm^ from the experiment [25]. 
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